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Abstract 

Motivated by the study of the potential use of blowing and venting operations of 

ballast tanks in manned submarines as an alternative control system for manoeuvring, 

we first propose a mathematical model for these operations. This model extends previ- 

,S^ ous works where only blowing is considered. Then, the model is applied to the control of 

rrt an emergency manoeuvre by using only blowing and venting. To this end, we formulate 

Cj a suitable constrained, nonlinear, optimal control problem where controls are linked to 

I I the variable aperture of blowing and venting valves of each of the tanks. The state 

law is composed of a system of nonlinear differential equations where the equations 

^^^ modelling blowing and venting processes are coupled with the Feldman, variable mass, 

f^ coefficient based hydrodynamic model for the equations of motion. In a second part, 

rsj we carry out a rigorous mathematical analysis of the model: existence of a solution 

\^ for both the state law and the optimal control problem is proved. Finally, we address 

^^ the numerical resolution of the optimal control problem by using a descent algorithm. 

f — Numerical experiments seem to indicate that, indeed, an appropriate use of blowing 

^^ and venting operations may help in the control of an emergency manoeuvre. 

^H 

^-H Keywords: Ballast tanks, manned submarines, blowing-venting operations, optimal con- 

^ trol, numerical simulations. 
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^ 1 Introduction 



Manned submarines are equipped with several ballast tanks distributed along its hull. When 
filled with water, they contribute with the submarine mass allowing it to submerge. During 
an unexpected event or emergence they act like a dispositive for emerging to the surface: air 
is blown into the ballast tanks from very high pressure bottles expelling the water out of the 
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tanks. The submarine loss weight, its buoyancy is higher, and it can emerge quicker. In the 
last years, several works have addressed these emergency rising manoeuvres (see pil^ fT^ITB] 
and the references given there). To fill the tanks with water again, air is vented out of the 
ballast tanks. A valve located at the top of each of the tanks is opened, air escapes outside, 
and water flows back into the tanks. 

On the other hand, under certain circumstances, like gathering intelligence missions or 
special operations, submarines may need to perform manoeuvres with very specific require- 
ments. In these cases, submarines often perform small blowing and venting operations, not 
because of an emergency, but to slightly modify the buoyancy of the vehicle. This way, 
blowing and venting become a complementary tool for manoeuvring. These manoeuvres 
are currently performed based exclusively on the operator experience and, due to the high 
degree of accuracy required, would enormously benefit from the implementation of a control 
system. 

The issue of modelling the blowing of ballast tanks has been addressed, for instance, 
in [2] and [13]. Up to the best knowledge of the authors, venting has not been addressed 
so far, and most importantly, the coupling of both processes as a control system has not 
been considered before. This is the first aim of the present work. In Section [2] we propose 
a model for a coupled system of blowing-venting operations. For the particular case where 
only blowing is considered, the model presented in this work was numerically compared in 
[7] with the one proposed by Watt [13] giving similar results, but also showing the ability 
to capture some phenomena that were overlooked by that last model. Then, these two 
processes are coupled with the usual Feldman's coefficient based hydrodynamic model for 
the equations of motion (see |B] and more precisely [IH]). We notice that, since the mass 
of the submarine changes with blowing- venting operations, some of the parameters (e.g., 
moments and products of inertia, weight, mass, etc..) that remain constant in Feldman's 
model are, in our case, time dependent. 

In a second part, we analyze the potential use of blowing- venting operations as a control 
system in a typical emergency manoeuvre. To this end, we model such a manoeuvre as a 
constrained, nonlinear, optimal control problem which has the aperture of blowing-venting 
valves of each of the tanks as the control variables. The underlying state law is composed 
of a nonlinear system of 3iV -|- 12 ordinary differential equations (ODE's), where N is the 
number of ballast tanks. 

The mathematical analysis of the model is presented in Section [4] and is divided into two 
parts. First, the existence of solution for the state law is obtained by using the classical 
existence theory for ODE's. However, due to the size and complexity of the mathematical 
model, checking the sufficient conditions for the existence of a solution is not an easy task. 
In particular, we carry out a detailed analysis of the main properties of the time- dependent 
mass matrix of the system, Subsection |4.1.1| Second, we also prove that the optimal control 



problem is well-posed, that is, that there exists a solution to it (Theorem 4.3 ). This existence 
result is obtained by using Filippov's theory for Bolza-type optimal control problems. The 
fact that controls appear in linear form simplifies the analysis. 

At this point, it is convenient to comment on the choice of an open- loop control system 
as an appropriate formulation for the problem considered in this work. Up to now, there 
is not an available technology to implement a control system as the one described in this 
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paper in a real manned submarine. In fact, the present study is motivated by the company 
Navantia S.A. Shipyard (Spain) as a first step to analyze the effectiveness of blowing- venting 
operations as a control system. In addition, the results provided by the formulation consid- 
ered here lets naval architects fit several parameters of blowing- venting systems (such as size 
of ballast tanks and/or blowing and venting valves) during the preliminary state of design 
of a prototype. Of course, once the capabilities of blowing- venting operations are positively 
evaluated by an open-loop control system, a second step would be the design of a closed-loop 
control system to be able to correct errors in the model, for instance in the estimation of 
the hydrodynamics coefficients, or the effects of disturbances acting on the surroundings of 
the vehicle. We plan to address this issue in a future work. 

Finally, for the numerical resolution of the optimal control problem we use a gradient 
descent method as in [TU]. The performance of the algorithm is illustrated through the 
numerical simulation of an emergency manoeuvre. Numerical simulation results seem to 
indicate that blowing- venting operations may help in a significant way to the control of this 
type of manoeuvre. 

To conclude it is important to point out that although the study of the present work has 
been only applied to a typical emergency manoeuvre, the scope of the model and the tech- 
niques developed here is not limited to this particular situation. Indeed, the approach pro- 
posed in this paper can easily be extended to analyze the use of blowing- venting operations 
as a depth controller in snorting condition, which is a relevant topic in naval engineering. 
The stabilization, based on blowing and venting of ballast tanks, of other marine systems, 
like offshore facilities or pontoon docks is another possible application. 

2 Mathematical modelling 

The three-dimensional equations of motion for an underwater vehicle are usually described 
by using two coordinate frames: the moving coordinate frame which is fixed to the vehicle 
and is called the body-fixed system, and the earth-fixed reference frame which is called the 
world system. The position and orientation of the vehicle are described in the world system 
while the linear and angular velocities are expressed in the body-fixed coordinate system. 
These quantities are defined according to SNAME notation [S] as 

V{t) = [vnt),rint)] , Viit) = [x{t),y{t),z{t)f, r^S) - m)Mt)Mt)f 
v{t) = [ul{t),vl{t)\ , u^it) = [uit), v(t),wit)f , V2(t) - W),q(t), r{t)f . 

Here t is the time variable, st^ stands for the transpose of vector a, r}i denotes the position 
of the vehicle in the world system, 7)2 is the orientation in the same reference system, Vi is 
the vector of linear velocities in the body-fixed frame (where as usual u is surge velocity, v 
is sway velocity, and w is heave velocity), and finally 1^2 is the vector of angular velocities 
in the body-fixed reference system (p is roll rate, q is pitch rate, and r is yaw rate). See 
Figure [l] 

The mathematical model for the equations of motion is based on Gertler and Hagen's [9] 
six degree of freedom (DOF) submarine equations of motion, which were revised by Feld- 
man [B]. Adapting these general equations to the particular characteristics of a prototype 
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Figure 1: Variables and coordinate systems. 



developed by the company Navantia S.A. Shipyard (Spain), a very similar coefficient based 
hydrodynamic model was analyzed in [10. This latter model will be the starting point for 
the more general model that we will introduce in this section. Anyway, to make the paper 
easier for readers we have collected these equations in Appendix [A} 



2.1 Blowing and venting model 

As shown in Figure ^ the blowing/venting system is composed of ballast tank, pressure 
bottle, blowing valve and venting valve. When the blowing valve is opened, air flows into 
the tank increasing the pressure and forcing the water to flow out through the flood port 
located at the bottom of the tank. When the venting valve is opened, air can flow out from 
the tank. The model can be divided into four main parts: 

• Air flow from pressure bottle. 

• Air flow through venting valve. 

• Water flow through flood port. 

• Evolution of pressure in ballast tank. 
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Figure 2: Blowing and venting processes. 



Finally, in Subsection |2.1.5| we shall introduce a new set of variables that will control the 
aperture of blowing and venting valves at each ballast tank, and shall summarize the whole 
system. 

Variables and symbols introduced in this section are summarized in Table [T] 



2.1.1 Air flow from pressure bottle 

When the blowing valve is opened, the air in the bottle is blown into the tank through 
a nozzle. Pressure losses and heat transfer in the tube that connects the bottle and the 
tank are, for the moment, neglected. However, as we will see later on, pressure losses can 
be indirectly taken into account. Under the above conditions, we need to study the one 
dimensional steady flow of an ideal compressible gas. This can be found in any classic text 
on fluid mechanics (see for example [1]) but for completeness we include it here. 

At the beginning of blowing, the flow is supersonic due to the high pressure difference 
between bottle and tank. As air flows out of the bottle, this difference decreases and the 
flow becomes subsonic at a certain time. This transition happens when pressures in bottle, 
Pp, and tank, ps, are such that 



= Pr 



where Pr 



^ 2 , 



with 7 the isentropic constant. 



Supersonic flow: In this case, it is easy to calculate the mass flow rate at the nozzle 
throat, where the Mach number is unity. Precisely, 



* * /I ra A 

'^F ^ Pa^ ^ = Pafiao A, 

Pa,0 ao 



(1) 
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A 


area in nozzle throat {m?) 


Ah 


outlet hole area [m?) 


A, 


vent pipe section {m?) 


Ch 


outlet hole coefficient 


Htk 


ballast tank height {m) 


hwcit) 


height of water column in the tank (to) 


msit) 


mass of air in baUast tank [kg) 


mpit) 


mass of air in pressure bottle {kg) 


mpit) 


mass flow rate from pressure bottle {kg/ s) 


rhy{t) 


mass flow rate through venting valve {kg/ s) 


PB{t) 


pressure in ballast tank {Pa) 


Pext{t) 


pressure outside the venting valve {Pa) 


PF{t) 


pressure in bottle {Pa) 


PSEAit) 


pressure outside the outlet hole {Pa) 


qsit) 


water flow through outlet hole {m? / s) 


Tb 


water temperature {K) 


Trit) 


temperature in pressure bottle {K) 


Vbo 


initial air volume in ballast tank {rm?) 


Vnit) 


volume of air in ballast tank {m?) 


Vbb 


ballast tank volume {m^) 


Vf 


pressure bottle volume (to^) 


{xb,yb,zb) 


location of geometrical center of ballast tank (to) 


Zh 


outlet hole distance from origin (to,) 


Zt 


tank top distance from origin (to,) 


Zy 


venting valve distance from origin (to) 


1 


isentropic constant 


P 


density of water {kg/m?) 



Table 1: Variables and symbols 

where A is the area in the nozzle throat, a is the speed of sound, v is the air velocity, pa is 
the density of air, the asterisk means conditions wherein the Mach number is equal to unity 
and the subscript denotes stagnant conditions. These two conditions are related by 



Pafl 
Pa 



7-1 



T* 



7- 



1 



where T is the temperature. 

Assuming that in this particular case stagnant conditions are the bottle conditions, using 
the above relations, the ideal gas law p^fl ~ J'°j, and ao = y^jRgTo, the mass flow rate is 



VR9Tp{t) V 



2 

7 + 1 



(2) 
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where Rg is the gas constant for air and Tp is the temperature in the bottle. 
Subsonic flow: In this case, the mass flow rate is given by 

mp = Pa.eMeUeAe = Pafltto-^^ — M^A^, 

Pa,0 iQ 



(3) 



where A,, is the exit area, in this case A^ ~ A, AI is the Mach number and the subscript e 
denotes conditions at the exit. The stagnant and static conditions are related by 



To _ f Pafi 
T ^\pa 



7-1 



= l + l^M\ 



Particularizing for po ~ pf and p = pe = pb , the Mach number at the exit is given by 



M, 



\7-l 



EL 

Pb 



(4) 



Substituting Q into ([3]) and proceeding as in the supersonic case, we get 

ApFit) 27 



mp(t) 



y/R^TW)\ 7 



PB{t)\-' fpB{t)\ -< 



PF{t)j \PF{t) 



(5) 



The mass flow rate from the bottle can then be expressed for both the subsonic and 
supersonic cases as 

(6) 



rnpit) = ^!'^^^\,, ■ p{pB{t),PF{t)) 



^RgTF{t) 
with p^{pB{t),PF{t)) = p,na. = Cte if ^ > ( ^) ^ . 

As the pressure in the bottle drops, the temperature increases. If the heat transmission 
is neglected, however, the process can be considered to be adiabatic. As mentioned in [2], 
experimental results sustain this assumption. Let mFQ^ Pfo be the initial mass and pressure 
in the bottle and Vp the bottle volume. Under the above assumptions the momentary 
pressure and temperature are given by 



PFit) 



Hence, 



mF{t) 
mpQ 



PF{t) 



PFO, Tpit) 



PF{t)VF 

RgiriFit) 



(7) 



rnlnVF 



\/RgTF{t) V '"-FO 

Thus the final equation for the mass flow rate is 

mFft)=Af "'^^y^'P^° ) p{pB{t),mF{t)) 
V m'pgVF 



(8) 
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with 



/ 2 \"'-i 



Pr< 



EL 

Pb 



n{pB,mF) = < 



\ 



27 
7-1 



PB 



PB 



,PFa 






,PF0 



\niFO J 



1<^<P. 



Pb 



0, 



EL 

PB 



< 1 



(9) 

Since the initial mass flow rate depends only on the initial conditions in the bottle, it 
can be considered as a constant with value 



771^(0) = A 7 



7 + 1 



-1+1 \ 



Vf J 



(10) 



This initial mass flow rate has been measured for several blowing intensities. Let rhFmax 
be this measured maximum mass flow rate. Instead of using the numerical value of area A 



when computing (|8|, we obtain it from (10) and the measured value for rriFmax, that is 

1 
Vf \ 



A — niFr, 



.7 + 1 






7P_F0"T.F0 



(11) 



By doing so, we ensure that the initial mass flow rate calculated coincides with the real 
measured value. This way, although pressure losses are not considered in the model, they 
are indirectly taken into account. 

2.1.2 Air flow through venting valve 

Equations for the air flow from the tank can be obtained analogously to the ones used for the 
mass flow from the bottle. In this case, the pressures to evaluate are the pressure in ballast 
tank, Pb, and the pressure outside the vent valve pipe system, p^xt = Patm + P9{z + z^ — 
XbSinO), where Patm is the atmospheric pressure, z is the vehicle depth, z„ is the vertical 
distance between the origin of the body-fixed frame and the venting valve, and the term 
Xb sin0 models the tank height variation with the vehicle pitch. Thus, (z + z^, — Xf, sin0) is 
the depth at which the exit of the venting system is located. Therefore, similarly to (l6]), we 
have 

AvPB{t) 



Jhy{t) ^ fl{pB, Pext) 



RgTBity 



(12) 



where Ay is the vent pipe section and 



^J' = fJ-n 



if 



Pext 
PB 



< 



Pext 
PB 



Although the geometric details of the venting system are not available and therefore 
pressure losses can not be directly calculated, an evaluation oi'p{pB, Pext) for several values 
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Figure 3: Curve fitting result for fi{lV) considering pressure losses. 



of the ^^^ ratio taking into account the different pressure losses along the system has been 
provided. In particular, the critical pressure ratio is found to be ( ^^^^^ J — 0.5184 and 



M„ 



0.3550. Let H 



PB 



- be the pressure ratio. Using a least squares fit of this data 



(see Figure [3]), the following expression for /i is obtained: 

P-max if n < Ucrit 



M(n) 



-66.97n2 - 52.70n + 119.70 

n3 - i7i.3on2 - 75.65n + 294.6O 



if n. 



<n< 1 



(13) 







if n > 1. 



The variation in the mass of air in the ballast tank is the difference between the mass 
flow rate from the bottle and the mass flow rate through the venting valve. This way, by 



(121 and (13), 



mB{t) = -liiFit) -Jl{U{t)) 



AvPBJt) 

RgTB[ty 



(14) 



2.1.3 Water flow through flood port 



The difference between the tank and outside pressure forces the water to flow in or out from 
the tank trough the flood port located at the bottom. Let us assume, in the first place, that 
the pressure in ballast tank is greater than the outside pressure. Then, the water flows out 
from the tank. A detailed analysis of a draining tank filled with an ideal fluid can be found 
in [4]. The pressure directly above the flood port is psit) + pgh^cit), where hwc{t) is the 
height of the water column in the tank. Assuming that the tank has a regular shape, the 
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height of the water column can be expressed as 



n-wcKt) = 1 ^ 1. tltk. 



''max 



where Htk is the bahast tank height, A?Ti(t) is the mass loss in the tank and l^mmax is the 
maximum value of the mass loss (see Section [2. 2| for more details). 
The Bernoulli equation applied at both sides of the port gives 

PB{t) + pghwcit) - Apioss = PSEA{t) + ^pvfiit), 

where Vh and Apioss are, respectively, the velocity and pressure losses in the outlet hole. Both 
major and minor losses have been calculated by a separate model and can be compressed 
into the one single term 

C,h is estimated to be C,h = 2.5. 

This way the velocity in the flood port is given by 



,,. h{PB{t) + pgh^c{t) -psEA{t)) 



P(l + 0.) 
and the volume flow from the ballast tank by 



,B{t) = c,A,.,[t) . aA. /(^^"(^^+^-';^7^;\'^^""^^» , (15) 

V py^ + '^h) 

where 

PSEA = Patm + Pg{z + Zh ~ Xb sin^). 

Here Zh is the outlet hole location. Ah is the outlet hole area and Ch is a coefficient that 
takes into account that, since the outlet hole is actually a grid, the effective area is smaller 
than Ah- 

If pressure in ballast tank is lower than the outside pressure, then water flows into the 



tank. In this case, the expression for the volume flow is (15), but with a negative value 



2.1.4 Evolution of pressure in ballast tank 

When the blowing valve is opened the air is blown into the tank at a very high velocity, 
rapidly mixing with water. This promotes good heat transfer from the water to the ex- 
panding air and we can work under the assumption that the air will immediately adopt 
the temperature in the tank. The process can then be considered to be isothermal. As 
mentioned in [2], experimental results sustain this assumption. 
From the ideal gas law it follows that 

. ,_,. mB{t)RgTB -VB{t)pB{t) 

^"^'^^ v^) ■ 



2 MATHEMATICAL MODELLING 11 



The volume occupied by air increases (and in the same quantity) as the volume occupied 
by water decreases. Thus, the rate at which it changes, Vg, will be the volume flow out of 



the ballast tank, qs, given by (15). Using the perfect gas equation, the momentary volume 



of air in the tank is Vb (t) = j-^ . The variation in tank pressure is then 

PB{t) 

. ,,^ PBJt) . , . P%it)<lB{t) , . 

PB[t) --mB{t) ^ . ,^ n ^ ■ (16) 

mB[t) mB{t)RgTB 

At the mathematical and numerical levels, the presence of the square root in the equation 
(15 1 generates serious difficulties. Indeed, if the term inside the square root vanishes, then 
the gradient blows-up. To overcome this difficulty we approximate the square root near the 
origin (x G [0,^]) by a fourth-order polynomial as follows. Consider the polynomial 

P {x) = ao + aix + a2X^ + a:),x^ + a4X . 

The coefficients of this polynomial are determined in such a way that the following conditions 

hold: 

(a) P(0)-0 

(c) P(0) = 

(e) Jo P {x) dx = j^ ^dx. 

Notice that conditions (a)-(d) are continuity conditions of P and its derivative P at the 
extremes a; = 0, ^, and condition (e) has the effect of preserving the average volume flow. 
Numerical tests showed that ^ = 1 gives a good approximation for the square root in 



( |T5| . We obtain 

flo = «! = 0, 02 = 8.75, a^ == —14 and 04 = 6.25. 
To sum up, we consider the new function 



P{x) 



~\J—x if a; < — ^ 

-P[-x) if -^<x<0 

P{x) if < a; < C 

yjx if X> ^ 



and replace qB (t) by its approximation q^ (t) defined as 

- f.. r< A p('^^PB{t)+ Pghn,c{t)-pSEA{t))\ ,^„. 

qB it) - C^A.P (^ -^^^^ j . (17) 

2.1.5 Blowing and venting controlled system 

Once a model for blowing an venting operations has been presented, our next goal is to use 
such a system as a control mechanism to improve the manoeuvrability of the vehicle. To 
this end, we introduce a new set of variables: the control variables of the blowing-venting 
system. Let Si, s^ G L°° {0,tf; [0, 1]) denote, respectively, the grade of aperture of blowing 
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12 





MBT 2 


MBT 3 


MBT 4 


MBT 5 


An 


0.191 


0.191 


0.191 


0.191 


A, 


0.0177 


0.0177 


0.0177 


0.0177 


Ch 


0.7 


0.7 


0.7 


0.7 


Htk 


5 


5 


6.2 


6.2 


PFa 


2.5 • 10^ 


2.5 • 10^ 


2.5 • 10^ 


2.5 • 10^ 


Tfo 


293 


293 


293 


293 


Vbo 


0.001 


0.001 


0.001 


0.001 


Vbb 


21.4 


21.4 


22.9 


22.9 


Vf 


0.8 


0.8 


0.8 


0.8 


Xb 


-28.6 


-28.6 


23 


23 


Vb 


1.2 


-1.2 


1.7 


-1.7 


Zb 


0.595 


0.595 


0.975 


0.975 


Zh 


-2.895 


-2.895 


-3.897 


-3.897 


Zt 


2.105 


2.105 


2.303 


2.303 


Zv 


2.705 


2.705 


2.705 


2.705 



Table 2: Ballast tanks characteristics. 



and venting valves of the i— th ballast tank during the time interval [0,i^]. Thus, from 



( 14 1 and ( 16 ) we obtain the final equations governing the controlled evolution of the mass 



of air in any pressure bottle {mp.{t)), the mass of air in the corresponding tank (7713. (i)), 
and its pressure {pBi{t)). We obtain 



rriFiit) = Si{t)A 



mF,(t)'y+V_Fo 



fJ'i{PBi{t),mFt{t)) 



(18) 



mB^{t) + mFtit) = -Aii(n(t)) ^— -^ : 

\/nglB 



(19) 



rriBiit) . 



i{t) ~mBi{t) 



PBiit)qBtit) 



(20) 



PBi(i)'^"""' ■ — V"/ RgTB 

The above formulation assumes that the flow trough the valves varies linearly with their 
aperture, Si(t), Si{t). Of course, once the control system were implemented in a real vehicle 
this assumption should be adapted to the particular characteristics of the chosen valves. 

The values of all the required geometrical parameters for the four tanks considered here 
are summarized in Table [21 
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2.2 Coupling of blowing- venting system with a variable mass model 
for the equations of motion 

As water flows in or out of tlie tanks there will be mass variations located at several points of 
the vehicle. Since the equations of motion assume the mass of the submarine to be constant, 
it is necessary to identify which terms, formerly constant, will become time dependent due 
to its dcpcndancc with mass. We will need to write the following properties as a function 
of the amount of water in the tanks: 

• Mass (to). 

• Weight (W). 

• Moments and products of inertia (/j,, ly, I^, Ixy, ■ ■ ■)■ 

• Location of the center of gravity (xq, j/g, zq)- 

Let us assume there are N ballast tanks with geometrical centers located at points {xbi, 
Vbi, Zhi) (where the subscript i denotes the «— th ballast tank). Let toq be the initial mass 
of the submarine (with all tanks completely filled with water) and Arrii the mass loss in the 
z— th tank. It is when the tank is completely filled with water, and reaches its maximum 
value when it empties. 

The volume of water that has left the tank is equal to the volume occupied by air except 
for the initial air volume in the tank, Vbo, which depends on the initial mass of air in 
the tank, mso, and the initial depth. The mass loss in the i—th tank can be obtained by 
multiplying this volume by the density of water, p, that is, 

Am,(t) =. p iVmit) Vbo) = P f "^^'iM^ - Vbo) . (21) 

V PBt[t) J 

The momentary mass of the submarine is then 

N 



TO(i) == Too - ^ Atoj (i) , (22) 

and the weight 



i=l 



W{t)^glmo-J2^^^m^ (23) 

with g the acceleration due to gravity. 

Although the vehicle buoyancy B is usually assumed to be constant, it is really a function 
of the vehicle depth. Indeed, as the depth increases, the outside pressure compresses the 
vehicle and therefore its volume decreases. Let Bq be the buoyancy at zero depth. The 
buoyancy is assumed to be linearly dependent with respect to the vehicle depth in the form 

.^..(.-^.) 

which means that the buoyancy decreases by a 1.5 percent a depth of 300 meters. 
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*z™, (t) 



//// 



Figure 4: Location of the mass loss. 



To find expressions for the location of the center of gravity and moments and products 

of inertia we will work under the assumption that the mass loss in each tank occurs at a 

point, (xbi, ybi, z„iii{t)), where Zmiiit) is the height at which mass loss happens for each 

tank. It varies, as shown in Figure l4J from the top of the tank, when it is completely filled, 

to its geometric center, Z{,i, when it is completely empty This variation is assumed to be 

linear so that 

{zti - Zbi)A.mi{t) 



Zmli ['' ) — Zf{ 



Am, 



where zu is the location of the tank top and Ami^max is the maximum value of the mass 
loss, i.e, the value when the tank is completely empty. 

Let 1^0, lyo, /^o, Ixyo, Ixzo, lyzo and xgo, 2/go, zgq be respectively the initial (all tanks 
completely filled) moments and products of inertia and coordinates of the center of gravity. 
Then the moments and products of inertia are given by 

N 

Ix{t) = {y^ + z^) dm = 4o - XI (Vbr + Zmuitf) Am,(t), 
•' i=i 

f ^ 

Iy{t) = {x^+ Z^) dm = IvO~Yl (^b^ + Zrnliit)^) Am,{t) , 

^ j=l 

N 

Iz{t) = J ix^+ y') dm = IzQ-Yl i^l + yl) ^"^^(i)- 

Z — 1 

N 

Ixy{t) = I xy dm^ Ixyo - ^ XbtybtAmi{t), 

i—l 

N 

Ixz{t) = xz dm = I^zo -y^^XbiZmH{t)Ami{t), 

2 — 1 

N 

Iyz{t) == yz dm = ly^Q ~^ybiZraii{t)Am,{t), 



(25) 



and the coordinates of the center of gravity by 
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rno - 


1 


mo - 


-Ef=iAm,(t) 
1 



[moycQ - Eili yu^miit) 



xcit) 

yait) 



zait) = -^i^ — 777 [moZGo - E^Ii Zjnh{t)Ami{t) 



"^0 - Ei=i Ami(^) 



(26) 



Substituting (22| — (261 into the equations of motion completes our variable mass model. 



2.3 Complete model in compact form 

Summarizing, the state variables of the system can be expressed in vector form as 

r "i"^ 

x(t) = \^mFXt).rnBXt),PBA't)\l<^<N ^'ni.tf' ,t^{if 

where N is the number of ballast tanks, r/ = [x, y, z, (f), 9, ip] is the vector of positions and 
Euler angles and v — [u, v, w,p, q, r] are linear and angular velocities. This means that the 
vector state x belongs to R'^^+^^. 

Moreover we write the control vector as 



u(i) = [si{t),Si{t)]^ 



<i<N 



Si and Si being, respectively, the aperture of blowing and venting valves of the i— th tank. 
Therefore, u G R^^. 

Finally, the state law is composed of equations ([l8l)-([20|) for the blowing- venting system, 
the six kinematic equations as given in Appendix |A.1[ and the six equations of the hydro- 
dynamic forces and moments (see Appendix A.2[ where the time variable parameters as 
described in Subsection 2.2 have been taken into account). In compact form we write these 
equations as 

A(x(t))x(i)-/(t,x(t),u(t)), (27) 

where A(x(i)) e 7\/((3Ar+i2)x(3N+i2) ^^^ /(t,x(t),u(i)) e RSA'+is for every i > 0. In 
Section [4] we will analyze in detail the structure of A and /. Nevertheless, at this point, it 
is convenient to comment on the explicit dependence of / with respect to the time variable 
t. Since we plan to analyze the potential use of blowing-venting as a control system for 
manoeuvrability, deflection of bow plane ((5(,(i)), deflection of stern plane {Ss{t)), deflection 
of rudder {Sr{t)) and propeller speed (n(i)), that typically are the elements used for the 
manoeuvrability of the submarine, will be fixed to some convenient values during the whole 
time of manoeuvre. More precisely, we assume that (S(,, 6s, Sr, n e L°°{Q,tf\K), where 



K = 



Svr 57r 
36' 36 



Stt 


57r' 


X 


36 


36 J 





111 Tvr 
36' 36 



[0, 2.5] 



Therefore, the explicit dependence of / on time is linked to these four functions (see Sub- 



sections 5.2 and A. 2 for more details on this passage) 
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3 Formulation of the control problem 

Next, we plan to analyze the potential use of blowing- venting operations as a control mech- 
anism in a typical emergency rising manoeuvre where the submarine must reach surface 
quickly while keeping its stability. We are now in a position to formulate the manoeuvra- 
bility control problem. Given an initial state x (0) = x° and a desired final target x*-f , the 
goal is to calculate the vector of control u = u{t), which is able to draw our system from 
the initial state x° to (or near to) the final one x*^ in a given time tf, also minimizing a 
cost functional. In mathematical terms we have the Bolza-type problem 

Minimize in u : J (u) = $ (x(i/),x*/) + / F{t,x{t))dt 

Jo 
subject to 

A(x(i))x(<) = /(<,x(i),u(t)) (Pt,) 

x(0) = xO e r^ 

< s,{t),s,{t) < 1, 1 <i< N, and x (i) e n, 
where ft stands for the set of constraints for the state variable. Typically 



3JV+12 

E 



$(x(t/),x*/)= Y. «.(^j(V)-^50 (28) 

with ttj > penalty parameters, and 



37V+12 

F{tMt))= E /3,(x,(i)-x,W)' (29) 

with f3j > also weight parameters and x(t) = [a;j(i)] a desired trajectory. 

The set fl, which models the constraints on the state variable, has the following structure. 
For the variables entering in the blowing-venting model, since the outside pressure at a 
certain depth is the sum of the atmospheric pressure and the weight of the water column 
above the submarine, even if the pressure in the ballast tank is slightly lower than the outside 
pressure and the vehicle is close to the surface, we can safely assume that the pressure in 
the tank will always be greater than the atmospheric pressure, that is ps. > Pg = Patm- 
Although an upper bound can not be so easily obtained, it is easy to see that the pressure in 
the tank will always take finite values, which justifies the assumption psi < Pg < +oo. It is 
also immediate to see that the upper bound for the mass of air in the tank and bottle is the 
initial mass of air in the bottle. By hypothesis, there will always be a residual amount of 
air in the tanks, mso- Therefore msi > Trig = mso > 0. As we stated before, the pressure 
in the tank will not drop below Patm- Since the air will flow due to the pressure difference 
between bottle and tank, the pressure in the bottle will have the same lower bound. This 
way, using the perfect gas equation and ([7]) a lower bound for the air mass in the bottle can 
be obtained, mp. > nip > 0. Summarizing, we are able to assume the constraints 

< rrip < niFi < rrip < +cx) 

< mj3 < niBi < m^ < +oo (30) 

< Pb ^ Pb, < Pb < +00- 
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As for Euler angles, since we are dealing with a manned submarine, typically 

TT , TT TT ^ Vr 

~-r<4><ir, --r<S<-, < ^ < 27r. 
4 4 4 4 

Due to the bounded nature of ocean, the position components (x, y, z) are also limited to 
some bounded rectangle. Finally, the physics of the problem also imposes a constraint on 
the rest of components (i.e, linear (u,v,w) and angular (p,q,r) velocities). 
To sum up, we can assume that i7 is a bounded rectangle. 



4 Mathematical analysis 

Once the mathematical model for the blowing-venting procedure has been established and 
the associated control problem has been formulated as a Bolza-type problem, then we will 
analyze it to provide an existence result. That is, we will show that, under suitable assump- 
tions, for any admissible initial state x" G il, there exists a control u(t) minimizing the cost 



functional in (P^ , ) 



4.1 Well-posedness of the state law 

We start our analysis by proving that for any initial state x*^ G fi and any admissible control 



u(i) there exists a unique solution of (27) starting from x*', and also satisfying x(i) G fi. 



defined in some interval < t <tf where tf = i/(x°) only depends on the initial condition. 
However, since the state law is expressed in implicit form, we can not apply standard results. 
To overcome this difficulty we will show that the matrix-valued map A is smooth and takes 
nonsingular values, that is, the state law can be rewritten in explicit form as 



x(i) = A(x(t))-i/(<,x(t),u(t)) 

4.1.1 A is nonsingular- valued 

It is clear that for any x G 17, the matrix A(x) has the form 



(31) 



A(x) 



BV(x) O3NX6 OsNxe 
OexsN h Ogxe 

Oex3N Ogxe M(x) 



(32) 



where BV(x) is the submatrix associated to the blowing-venting equations and M(x) is the 
so called (variable) inertia matrix. Furthermore, BV(x) is a 3iV x 3iV matrix structured in 
3x3 diagonal blocks 



BV(x) = 



BVi(x) 







3x3 



J3x3 



BVjv(x) 



where, for 1 < i < A^, 



BV,(x) = 



1 








1 


1 








-1 


niB 



PBi 



(33) 



(34) 
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withdetBV,(x) 



> —f- > by ( 30 ) . Thus B V^ (x) is nonsingular with 



BV,(x) 



1 








-1 


1 





-PB, 

rriB, 


TIB,- 


PB, 
TUB 



(35) 



and, therefore, the fuU matrix BV(x) is also nonsingular for every admissible state x G f2 

with 

BVi(x)-i ••• 0, 



BV(x)-i 



-"3x3 



(36) 



03x3 ••• BVAr(x)-l 

It just remains to prove that the inertia matrix M is invertible. From the dynamic 



equations of motion and (22)-(26l we have that M(x) has the form 

M(x) =M„(x)+Mc 
where 



M,(x) = 
is the variable part of the matrix, with 





■m{x)l3 -S{x) 
5(x) /(x) 



(37) 
(38) 



5(x) = 



-m(x)zG(x) m(x)yG(x) 
m{x.)zQ{x) — m(x)a;G(x) 

~to(x)2/g(x) to(x)xg(x) 



and 



Ixix) -Ixyi^) -Izxi^) 

7(x) = -4y(x) lyix) ~Iyz{^) 

-Izxi'X.) -Iyz{^) Izi'X.) 

the inertia tensor, and Mc is the so-called added inertia matrix 



PlSv/ 
2' ^u 




















2' ^v 





2' -^p 





P;4v' 
2' ^r- 








-'^Pz'^ 





p/47/ 
-2* ^g 








2' ^^ii 





2' "p 





-fz^i^; 








2' -"^li 





-¥^M'. 








2' ^^i] 





2' ^^p 





-¥^K 



(39) 



It is usual in the literature on dynamics of submerged vehicles (see.S , Property 2.4, for 
instance) to assume that M^ is a symmetric and positive definite matrix (and therefore 
invertible). However, experimental values of the non-dimensional hydrodynamic coefhcients 
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reported by Navantia showed that this is not a reahstic assumption in all cases. In partic- 
ular the inertia matrix used in our numerical experiments, based in the experimental data 
provided by Navantia, is not symmetric, but it is invertible (see Table p^ in the Appendix). 
Let Mq be the rigid-body inertia matrix of the submarine with all the ballast tanks 
completely filled with water, i.e., 



M 



— 



mo /a 



mo 





-zgo 
Vgo 



zgo 


-xqo 



-yco 

xgo 




ma 





-yco 



-ZGO 



XGO 



VGO 

-xgo 




(40) 



This matrix is symmetric and usually is assumed to be positive definite ([5], Property 2.2). 
Therefore whenever matrices Mc and Mo are assumed to be symmetric and positive definite, 
the inertia matrix Mo -I- Mc is invertible, since it is also symmetric and positive definite. 
Otherwise the invertibility of Mo + Mc must be checked for any particular model. 
Finally, nonsingularity of M(x) follows from the next Lemma. 

Lemma 4.1 Let D e TVf™^™ be an invertible matrix and let A e ^'"><™ be a square matrix 
such that 

|lD-iA||=sup{|D-iAy|:|y| = l}<l 

where I • I denotes the Euclidean norm in R™ . Then D + A is invertible with 



(D 



^(-l)«(D-iA)"D-i 



n=0 



Computing the operator norm of a matrix is not in general an easy task, but next Lemma 
provides a very useful estimate. 



Lemma 4.2 Let D = [dij] e M' 



be a square matrix. Then 

I|d|I<,/E4- 



(41) 



The right-hand side in (41) is called the Frobcnius norm o/D, ||D||i?. 



Let X e f2 be an admissible state. Since M(x) = (Mo + Mc) -I- (M„(x) - Mo) with 
Mo + Mc invertible, by Lemma |4.H it suffices to show that 



(Mo + Mc)~'(M,(x)-Mo) 
to have the nonsingularity of M(x). Furthermore, Lemma 



< 1 

also ensures that 



4.1 



M(x)-i = ^(-1)" ((Mo + Mc)-^ (M„(x) - Mo))" (Mo + Mc)"^ 



(42) 



(43) 



n=0 
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Summarizing, M(x) is invertible, for any state x G ft, whenever Mg + Mc is invertible 
and the variation of mass due to the blowing-venting manoeuvres is small enough. More 
precisely, if 

m^ < —^^^^^^^^=^^^^^^^^= (44) 

holds, where Am^ = X]i=i ^"^i.maa: — P^i^ii^BBi — Vbo) is the maximum variation of 
mass, Xmin > is the smallest singular value of the inertia matrix Mq + Mc and c^ is the 
geometrical center of the i— th ballast tank, then A is nonsingular- valued. 



Theorem 4.1 If the matrix NIq + Nlc is invertible and the inequality {44') is satisfied, then 



the matrix-valued map A : Q — > _y\/((i2+3Af)x(i2-i-3Af) i^f^gg nonsingular values, that is, A(x) 
is invertible for any x G fi with 



A(x)-i = 



BV(x)-i Oawxe Osnxq 
Ogxsn I& Ogxe 

06X3W 06x6 M(x)-1 



Moreover, the map x ^> A(x) ^ is continuously differentiable. 



Proof. To show that A(x) is invertible for any x e J^, it suffices to show that (42) is satisfied 
uniformly with respect to the state variable. We have 



(Mo + Me)-i(M,(x)-Mo 



< 



An 



(Mo + M,)-i ||M,(x)-Mo| 
^ '|M,(x)-Mo|| 



(45) 



by properties of the matrix norm (see Appendix A in |llj . for instance). Let us now obtain 



estimates for the supremum of the Frobenius norm of the variable matrix. Firstly, from ( 21 ) 

N 



^Ami(x) 



N 



N 



< ^ Anii^jnax = P^iVBBi ^ Vbo)- 



Let us now estimate the variation of the coordinates of the center of gravity 



N 



^Xh,Am,{x.) 

N 

^ytiAmi(x) 



N 



< (Am+)^|xb,|, 

N 

<{Am+)Y,\yM\ 



4=1 



and 



N 



^z„H(x)Amj(x) 



N 



< 



(Am+) J2 



Zbi 



(46) 

(47) 

(48) 
(49) 



with {xbi,ybi, Zbi) the coordinates of the geometrical center of the i— th tank (see Figure [4]). 
Finally, from ptl) and pel) 



|4(x) - ho\ = 



N 



J2 iyb^ + ^nih) ^rn,{yi) 



N 



< (Am+) Y: {yl + zl) 



(50) 
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and, in a similar way, 

|/,(x)~/,o|<(Am+)^ 
and 



N 



i=l 



N 



l^bi + ^bi) 



|4(x) - U < (Am+) J2 (4 + yl) (51) 



4 = 1 



AT 



AT 



|/,,(x)|<(Am+)^ 



Xbi\ \yu\ 



\Ixz{'x.)\ < (Am+) ^ |a;f,i||zhj 



AT 



1=1 



(52) 
(53) 



Hence, from Lemma 4.2 and estimates (46l-(53), it follows that, for any x G f7. 



lM,(x)-Mo|l<(Am+) 



\ 



i + mlf24^ + yl + 4^ 



(T.^l + yl + '^lj 



(54) 



where we have also used the convexity inequality ( X]i=i '^0 ^ ^ [ X]i=i '^1 ) ^^'^ ^^^^ t^^^t 
a5 < (a^ + b'^)/2. Combining this inequality with (44), we have the desired inequality (42). 



Finally, from (36) and (43) it follows that x ^> A(x) ^ is continuosly diffcrcntiable and the 
proof is complete. 

Remark 1 For the particular data of the prototype considered in our numerical experiments, 
Mq + Mc is invertible. Moreover, Xmin = 1-54 • 10^^ and hence 



1.73 -10^^ (55) 



(Am+) J 3 + AN (Zl, xl + yl + zl) + | (^l, x^ + yl + zl) 



4.1.2 The state law is well-posed 

Let K, = [0, 1]^^ be. An admissible control of the problem will be an essentially bounded 
map taking values in /C, that is, L°°(0, +cx);/C) is the space of admissible controls. Given 



u g L°°{0,+oo;IC), our aim in this section is to show that the system (31) has a solution 
for any initial state x°. We begin by recalling the classical theory on this subject. For the 
details we refer the reader to [11] Appendix C]. 

By a (Caratheodory) solution we mean an absolutely continuous function t ^> x(^) G Q, 
defined on some interval / = [0,t/], which satisfies the integral equation 

X (t) = x*' + / g{s, x(s)) ds, for every t € [0, tf] , 

"'0 

where we write g : I x il — > R'^^+-'^^, g{t,x) = A(x)^^/(i,x, u(i)) for simplicity. 

As it is well-known, if g satisfies conditions (H1)-(H4) below, then we can ensure the 



existence and uniqueness of a maximal solution of (31 ) for any initial state: 
(HI) for each x G J7, the function g (•, x) : / — > R^ is measurable. 
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(H2) for each t & I, the function g (t,-) : Q ^ R^ is continuous, 

(H3) g is locally Lipschitz with respect to x, that is, for each x" € fi there exist a real 
number e > and a locally integrable function a : I ^ R+ such that the ball B^ (x°) 
of radius e centered at x'^ is contained in 51 and 

\g{t,x)-g{t,y)\<a{t)\x-y\ for every i e / and x,y eS^ (x") , (56) 

(H4) g is locally integrable with respect to t, that is, for each x" e fi there exist a locally 
integrable function /3 : I ^^ R+ such that 

|g(t,x°)| </3(i) a.e. t e I. (57) 

Our main result in this section follows: 
Theorem 4.2 For each x*' e 51 and each u G L°° (R"'";/C) there exists a unique maximal 



solution of (31) defined in [0,t/] , where tf = t/(x'^) only depends on the initial state x'' and 



is uniform with respect to the control u £ L°° (R+;/C) . 

Proof. To begin with, we notice that thanks to Theorem |4 . 1 1 we may restrict the analysis 



that follows to the right-hand side of the state law (27). Let us see that conditions (H1)-(H4) 
above hold. Since the time variable t only appears in the control functions (s^ (t) ,Si (i)) , 
1 < i < N, and, by hypothesis, these functions belong to L°° (R"*"; [0, 1]) , it is clear that for 
each X e f2, the function t "^ f (t, x) is measurable. 

As regards the continuity of the function / (t, •), only the equations modelling the air 
flow from pressure bottles need of a detailed analysis to check that the pass from supersonic 
flow to subsonic one is continuous. A direct computation shows that this is so whenever 
mpi > rup > and ps^ > Pg > 0- 

As for condition (H4) , from the form in which controls appear in the state law it follows 
that for each x° S 51 the components fi (t, x") , 1 < i < 3A'' + 12, of / (i, x°) are uniformly 
bounded with respect to i, that is, ( [57| ) holds for a constant function /3 (t) = (3 and what is 
more important, this constant is uniform with respect to u € L°° (R+; K.) . 

Next we analyze the local Lipschitz condition of / (i,x). Consider firstly the equations 
/i, i = 1 + 3j, < j < iV — 1, which model air flow from bottles. Taking into account the 



set of constraints (30), a direct computation shows that fi (i, •) S W^'^ (51) and therefore 
they are Lipschitz. As before it is important to notice that the estimates on the partial 
derivatives of fi {t, •) are uniform with respect to i, that is, there exists i > such that 



^f^ ,, . 



< L for every x G H and uniformly w.r.t. t > 0. (58) 

Similarly, functions fi,i = 2 + 3j,0<j<N~l, which appear in the equations modelling 



air flow through venting valve are continuous and satisfy an estimate as in (58). Thus, 
they are Lipschitz. As before, the Lipschitz constant is uniform with respect to the control 
variable. 

Consider now functions fs+sj, < j < N—1, which appear in the equations for evolution 
of pressure in ballast tanks. In this case, when pressure in a ballast tank equals outside 
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pressure, the derivatives of velocities in the corresponding flood port blows-up because of 
the presence of the square root. To avoid this singularity we have approximated the square 
root as shown in Subsection |2.1.4[ As a conclusion, once again we obtain that these maps 
are Lipschitz. 

The components fi, 37V+1 < i < 3A'^+6, only include the transformation matrix between 
body and world reference frames. Taking into account the constraints on Euler angles, it is 
clear that fi e C°° (fl) , 3iV + 1 < * < 3-/V + 6, and therefore they are Lipschitz. Notice that 
the time variable does not appear in these functions. 

As for the remaining fi, 3N + 7 < i < 3N + 12, these components include: (a) polyno- 
mial terms and terms in the form of absolute values; all of them are Lipschitz, (b) terms 
like Xj^/x'^ + x^. and \xj\ wx^ + x^ for some 18 < j < 24. Since these functions are contin- 
uous and the discontinuities of its derivatives are of a finite jump, they are also Lipschiz. 
(c) Since our model is mass variable, it is necessary to look carefully at the terms includ- 
ing mass m, weight W, center of gravity (xctVczg) and moments and product of inertia 
Ix,Iy, Iz^Ixy, ■ ■ ■ , because now they depend on some components of the state variable. Tak- 



ing into account the constraints (30) it is clear that these components are also Lipschitz. 
Finally, the product of Lipschitz functions is also Lipschitz. As before, these components 
fi do not include control variable u (i) and therefore the corresponding Lipschitz constants 
are independent of u (t) . 

This analysis lets us conclude that for each x" e fi and u eL°° (R+;/C) there exists 
a maximal time tf = t/(x'',u) and a unique maximal solution defined in [0,tf (x°,u)] . 
Looking at the proof of this existence result (see [11, Th. 36, pp. 347-351]) we realize that 
tf depends on u through the functions a{t) = a (u (t)) and f3 (t) — (5 (u (t)) which appear in 



(56 1 and ([57|. However, as shown before these functions a and /3 can be chosen uniformly 
w.r.t. u G L°° (R+;/C). As a conclusion, tf = tf (x°) only depends on the initial condition. 
This completes the proof. 

4.2 Existence of a solution for the optimal control problem (Ptf) 
Fix x'^ E n and let tf = tf (x*^) be the maximal time, as given by Theorem 



4.2 



for which 
system (|31[) is well-posed. The goal of this section is to prove the following existence result: 



Theorem 4.3 Let x'^ E fl and let < tf — tf (x'^j < +oo be as above. Then, there exists. 



at least, one solution to (Ptj )■ 



Proof. The proof follows as a consequence of the classical Filippov existence theorem 
for Bolza-type optimal control problems (see [H Th. 9.3.i, p. 314]). Indeed, let us see that 
the sufficient conditions for the existence of a solution hold. Due to the constrains on the 
state variable, the set A = [0,t/] x il is compact. Similarly, since the set /C = [0, 1] for 
control constraints is compact, the set A x JC is also compact. In addition, the function 



$ : ^ X ^ — > R defined by (28) is continuous. Continuity also holds for the functions F and 
g both defined on Ax IC. 



Notice also that, by Theorem 4.2 the set of admissible solutions for {Pt^ ) is not empty. 



Finally, we must check that for each (t, x) S .4 the orientator field 

Q {t, x) = { {z°; z) E Ri+3^+12 ■.z°>F{t,x), z^g {t, x, u) , with vl E K.) 
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is convex. This convexity easily follows from the facts that the control variable u appears 
in a linear form in the state law and the set /C is convex. 



5 Numerical analysis 

5.1 Algorithm of minimization 



There are several optimization methods which can be applied to solve {Pt^ ). Due to the 
complexity of the state law and the large number of variables involved in the problem, it is 
quite reasonable to use a gradient descent method with projection. Briefly, the scheme of 
this method consists of the following main steps: 

1. Initialization of the control input u*^. 

2. For fc > 0, iteration until convergence (e.g. |j(u'^+^) — j(u'^)| < e|j(u°)| , with 
e > a suitable tolerance) as follows: 

2.1 we consider the vector 



v^+i = u^ - AV J (u^) 
where A > is a fixed step parameter, and V J (u*^) is the gradient of the cost function. 

2.2 Since v*^"*"^ may be not admissible, we compute its orthogonal projection onto the 
admissibility set K, the unit rectangle in R*, that is, 

u'-^+i = Vk (v^^'+i) 

where 



u*^+i = min 



(max(0,v)-+i),l). 



The crucial step is the computation of the gradient V J (u*^) . This can be obtained by 
using the adjoint method which is described next: 

1. Given the control u^, k > 0, solve the state equation 

A(x(t))x(t) = /(f,x(i),u'=(t)) 
x(0) =x'=(0) 

to obtain the state x'^^^ (t) . 

2. With the pair (u'^' (t) ,x'^^^ (t)) , solve the linear backward equation for the adjoint 
state p (i) 

A (x^-+i it)fp (t) = -VxF (i, x'^+i (t) , u'= (t)) - [v./ (t, x^-+i (t) , u'^ {t))f p (t) 
A (x^-+i itf)fp [tf) = V.$ (x'^+i [tf) , x*/) 

where Vx is the gradient with respect to the state variable x. Thus, we obtain p'^^^ {t) . 
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3. Finally, 

V J (u'^) - V„F (i, x'^+i (i) , u^- (t)) + [Vu/ (t, x^-+i (t) , u'^ {t))f p'^+i (t) 
where now Vu is the gradient with respect to u. 
Here A"^ stands for the transpose of A. We refer to [5| for more details on this method. 

5.2 A numerical experiment 

In order to test the proposed models, this section shows the results of the numerical sim- 
ulation of an emergency rising manoeuvre. These results are used to analyze both the 
mathematical properties and the possible real-world applications of the proposed scheme. 
At each iteration of the gradient algorithm, the numerical resolutions of the state and ad- 
joint state equations have been carried out by using the ODE45 Matlab function, which is 
a one-step solver based on an explicit Runge-Kutta method. 

On an emergency situation, like on board fire or flood, submarines may need to rise 
to the surface quickly. The typical protocol for this situation is to use the control planes 
to pitch the nose of the submarine up, increase the speed, and blow the ballast tanks to 
reduce the weight of the submarine and drive it to the surface with buoyancy. As analyzed 
in [TJ I13j , small to medium size submarines exhibit a roll instability during these kind of 
manoeuvres. Full scale trials showed that submarines rose with roll angles up to 25 degrees. 
It is also known that if the the submarine emerges with a high roll angle, it may experiment 
large roll oscillations on the surface. Of course, this is an undesirable situation, particularly 
if the operators are attending to the original problem that required the emergency rise. 

For these reasons, it might be interesting to check if this situation can be prevented by 
using our control algorithm. To this end, an emergency rising manoeuvre has been simulated 
for two different scenarios: 

1) A standard manoeuvre. Initial depth is 100 m, initial speed is 2 m/s. At i = the 
stern and bow planes are set to —20 and 20 degrees respectively and the propeller is 
set to 150 rpm. At the same time, ballast tanks 2-5 are simultaneously blown with 
half the maximum intensity (this corresponds to s; = 0.5 in our model). Vent valves 
remain closed (s^ = 0) throughout all the simulation. Simulation ends when submarine 
reaches a depth of 10 m (an arbitrarily low value for which we can assume that the 
vehicle has reached the surface). To sum up: 



yt 



We note that, although we refer to this scenario as standard manoeuvre, the constant 
value for the deflection of control planes is of course a simplification of what would be 
done in real operation. 



Sr{t) = 







Ss{t) = 


= -20 




Shit) = 


20 




n{t) = 


= 2.5 




Siit) = 


= 0.5, i = l. 


..4 


Mt) = 


0, i = l. 


..4 
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2) Same manoeuvre, with the control algorithm acting from i = to t = 30 s using 
Scenario 1 as initialization and looking to achieve three main objectives: 

• Submarine must rise in a similar time as it does in the standard manoeuvre. 

• Rising pitch angle must be around 20 degrees and never above 25 degrees. 

• Roll angle must be as close as possible to zero throughout all the simulation. 

To this end, the following set of parameters is used: 

x(0) = {[mFo.,mBOuPBo^i<,<i, 0, 0, 100, 0, 0, 0, 2, 0, 0, 0, 0, o) 

t/ = 30 s 

z*/ = x\i, = 75, 61*/ = x*/7 = 10. 

0(t) = xi6(i)=O Vte[0,30]. 

ai5 = air = 1, "j =0, j ^ 15, 17 

/3i6 = 5, /3, = 0, jV 16 

A = 0.01 e = 10-4 

where the initial mass of air in the bottles is mpOi = 237.8376 kg, the initial mass of 
air in the tanks is msoi ^ 0.0126 kg and the initial pressure in the tanks is psOi = 
1.0846 • 10^ Pa. After t = 30 s, simulation continues with fixed controls (values equal 
to those used in Scenario 1) until the submarine reaches a depth of 10 m. 

Results for Scenario 1 (dashed lines) and Scenario 2 (solid lines) are shown in figures [5] to 
|9] Figure [5] shows vehicle depth, pitch angle and roll angle for both scenarios. For Scenario 2, 
pressure, mass of air in the tank and mass of air in the bottle are plotted for each of the tanks 
on Figure [6j The rest of state variables have not been included since they are not directly 
relevant for this particular manoeuvre. Comparison between dashed (standard manoeuvre) 
and solid (optimal controls) lines shows that the three objectives have been achieved for 
Scenario 2: the rising time is only a few seconds greater than in the standard manoeuvre, 
the final pitch angle is close to 20 degrees and the roll angle has been significantly reduced 
with respect to the standard manoeuvre. Indeed, Scenario 1 exhibits roll angles in the range 
of 3-4 degrees during most of the simulation time while in Scenario 2, after an initial peak 
of 2 degrees, the roll angle is kept within extremely low values during the whole t — [0, 30] 
interval. 

The optimal (solid lines) and standard (dashed lines) controls are shown in Figure It] 
(blowing valves) and Figure [s] (venting valves) . Although the convenience of the use of 
venting during an emergency rise can be arguable and is indeed not an usual practice, it is 
used here to demonstrate the algorithm capabilities. Very similar, although slightly worse 
results, were obtained by using only the blowing valves as control. 

As we can see in Figure [Tj the rolling moment is compensated by blowing more ballast 
from the starboard tanks. It may surprise that tanks 2 and 4 are being vented while there 
is no blowing air yet, but this seems to allow a smoother transition when the two valves 
are simultaneously open. As we said before, the use of venting valves seems to improve the 
results obtained by only blowing the tanks. 

The value of the cost function at each iteration is plotted in Figure |9] As we see, the 
algorithm shows exponential convergence. As is usual in this type of algorithm, results 
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Figure 5: Depth, pitch and roll angle for scenarios 1 (dashed line) and 2 (solid line). 
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Figure 6: Pressure (solid line), air mass in bottle (dashed line) and air mass in tank (dash-dot 
line) for tanks 2—5. 





Figure 7: Blowing valve aperture for standard (dashed lines) and optimal (solid lines) con- 
trols. 
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Figure 8: Vent valve aperture for standard (dashed lines) and optimal (solid lines) controls. 
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Figure 9: Cost function at each iteration. 
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depend on the initialization. Different initializations were tested obtaining different optimal 
controls. This seems to show the existence of several local and/or global minima. 

A Appendix: equations of motion 
A.l Kinematic equations 

(p ~ p + r cos((/)) tan(6') + q sm{(p) tan(6') 
9 ~ qcos((j)) — rsin((/)) 
r cos(0) + q sin((/)) 
^^ cos(6l) 

X — u cos(0) cos(V') + v(sin{(f)) sin(0) cos('0) — cos{(p) sin('0)) 

+ w(sin(0) sm{tp) + cos(0) sin(6') cos(V')) 
y ~ u cos{6) sin(-0) + v{cos(4>) cos{ip) + sin((/)) sin(0) sin(?/')) 

+ w{cos{(j)) sin{6) sin(?/;) — sm{(f)) cos(f/')) 
z — —u sin(6') + V cos{6) sin(0) + w cos{0) cos((/)) 

A. 2 Dynamic equations 

AXIAL FORCE EQUATION: 

m[u — vr + wq - xc{q'^ + r'^) + Vcipq — r) + zc{pr + q)] 

+ ^l^[x's^,y5l + x's^suHl + X's^sy^l] -{W-B) sin(0) + pT{l - tp) 
LATERAL FORCE EQUATION: 

m[v — wp + ur — yair^ + P^) + zciqr — p) + xc{qp + r)] 

- I' 



-I'^iY^r + Ypp + i;|,,|r|r| + Yp^pq] 



+ ^l^[Y^ur + Ypup + Yiv + Y^^wp] 

+ ^f[YiuHr + YiySrir, - ^)C] 
+ f ^'n'»w«w + {W-B) cos(0) sin(0) 
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NORMAL FORCE EQUATION: 



m[w - uq + vp — zaip +q ) + xcirp - q) + yc{rq + p)] 

.21 
2' i-qi ■ -qlqlilil ■ "rr- J ■ 2" 



= 7^1'^iK'i + Kkl'i^l^ + ^'r^^ + n^^[^'w^ + Z'^uq + Z'^pvp + Z'^^vr] 



Pl2, 



+!^r'[zy + ziuw + z',y] + ^p[z\^\u\w\ + z;^jv \M («' + «^')^ 



+i'' 



Zg^u 5s + -^a."" ^b + -^5,»7" (^s h? - 



C 



+{W - B) cos{e) cosicj)) 
ROLLING MOMENT EQUATION: 

IxP + {Iz - Iy)qr - hxT - IzxPq + lyzr^ - lyzQ^ + IxyPr - Ixy^l 



myow — myauq + myGvp — mzav + mzawp — mzcur 



Pl5 



P,5 



P,5 



Pl5, 



Pl5, 



j^Kpp+Ll^K.f+^l'K^^qr+'-l'K^^^^plpl + ^-l^K,^, 
+ ^l%up + ^l^Klur + P-1'k',v + ^-I'Kl^wp 



PjSj 



^-PKluv+^fK',,,,v\v\ 



^I'^K^u -r -(, i^^uu -r-i jv^i^i^,^, , 2 



^,3, 



,M- + -rj\„uv + -r A,,u,iu|u| + -i^K^^Sr - pQ 

+ {YgW - YbB) cos{9) cos(</.) - {ZgW - ZgE) cos{9) sin(</.) 

PITCHING MOMENT EQUATION: 

lyq + {Ix - Iz)rp - {p + qr)Ixy + {p^ - r^)Izx + iqp - r)Iyz 
+m[zG{u — ur + wq) — xg{w — uq + vp)] 



P,5 



P,4r 



P[M^q + M,.prp + M^^^^q\q\ + M,y] + 9 [M^w + M^uq + M^^^vr 



+ 'y 



My + M^uw + M^y + M^i^i^z 



y + u;^)5 



P,3 



r[M^^vw + M|„|u|u;| + Mywy + wy\ 



+'-p 



Ms u'5s + M^ySb + Ms yss V - 



c 



-{xgW - xbB) cos(6I) cos((/.) - (zgW - zbB) sin(6') 
YAWING MOMENT EQUATION: 

hr + {ly - Ix)pq - {q + rp)Iyz + {q^ - p^)Ixy + (rq - p)Izx 
+m[xG{v — wp + ur) — yG{u — vr + wq)] 
= ^l''[N'rr + Kyrl + Npp + N^^pq] + ^/^[iV^up + iV>r + iV>] 
Pii 



+y[Ny + Kuv + N^^^^j^v y + wy j 



+'y 



Nsy5r+Ns^y5rU--^\C 



P,3 



l-'N 



^^j^VW 



+ ixGW - xbB) cos{9) sin((/)) + (ycVK - ysS) sin(6') 

The value of the hydrodynamic coefRcicnts present in the above equations have been 
obtained experinientahy by using a scale model. 
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/ 


67 m 


Too 


2352 • 10^ kg 


Bo 


23073 • 10^ N 


Xgo 


Om 


Ygo 


m 


Zgo 


-0.264 m 


Xb 


m 


Yb 


m 


Zb 


-0.595 m 


IxO 


16390 • 10^ kgm^ 


lyO 


659770 • 10^ kgm^ 


ho 


659770 • 103 kgm^ 


IxyO 


kgm^ 


IxzO 


Okgm^ 


lyzO 


kgm^ 



Table 3: Submarine param.eters 



Kto 


0.5382 


Ktj 


-0.3811 


Kt.J2 


-0.1645 


KtJ3 





Kt.J4 





V 


1.09431 


c 


0.8976 


Kqo 


0.0701 


Kqj 


-0.0273 


KqJ2 


-0.0345 


Kqj3 





Kqm 





tp 


0.148 


Wf 


0.344 


D 


3.821 m 







Tabic 4: Propeller model eoefficients 
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K 


-0.46 • 10-3 


n 


-0.3 - 10-3 


^4 


-0.07- 10-3 


^rp 


0.6-10-3 


n 


-0.83 - 10-3 


z'^ 


-14.40 - 10-3 


X' 


1.42 • 10-3 


Y! 


-18.37-10-3 


Zyp 


-18.37- 10-3 


K\q\ 





Y' 
p 


-3.05 - 10-3 


Z'q 


-6.99-10-3 


^rr 


2.08 • 10-3 


Y' 
pq 


-0.07-10-3 


Kh 


-3.98-10-3 


^vr 


22.37-10-3 


Y' 


14.40 - 10-3 


z„ 


-3.96-10-3 


^uu 


-1.124- 10-3 


Y' 


7- 10-3 


Zvr 


-45.13 - 10-3 


^vv 


17.46 - 10-3 


Y' 


-4.65 - 10-3 


7' 


0.185 


VI 

WW 


7.75 - 10-3 


K 


-61.36-10-3 


Z'w 


-20.28 - 10-3 


-^w\w\ 





n%i 


-28.14- 10-3 


Z\w\ 


-1.99-10-3 


^wq 


-13.16- 10-3 


yi 


-0.83 - 10-3 


7' 

WW 


20.00 - 10-3 


^'w\q\ 





Y! 


0.67-10-3 


Zi 


-5.12- 10-3 


xu 


-3.90 - 10-3 


y: 





Zs,v 


-0.45 - 10-3 


xu 


-1.19- 10-3 


YvwN 







-5.12- 10-3 


^ks. 


-2.99 - 10-3 


Y' 


-0.1162 


z: 


-0.3-10-3 


^ks. 


-2.99 - 10-3 


Y' 


-0.1162 


z: 


-0.3-10-3 


M'. 


-0.98 - 10-3 


% 


-0.03-10-3 


7' 

^w\w\R 





M4 


-1.39- 10-3 


N( 


-1.20- 10-3 


Kp 


-0.07-10-3 


m;^ 


1.12 - 10-3 


K 


0.68 - 10-3 


K', 


-0.05 - 10-3 


m; 


-3.89 - 10-3 


K 


-0.68 - 10-3 


K 


-0.60- 10-3 


Kh 


-3.03 - 10-3 


Kq 


-0.91 - 10-3 


K 


-0.62 - 10-3 


M'„ 


-1.19-10-3 


K 


-4.83-10-3 


K' 

^p\p\ 


-0.30-10-3 


Mir 


-15.73- 10-3 


K\r\ 


2.09 - 10-3 


Kp 


0.30-10-3 


M^„ 


34.61 - 10-3 


K 


-18.64- 10-3 


Kr 


-0.21 - 10-3 


m; 


4.78 - 10-3 


N' 


18.50-10-3 


K 


0.26-10-3 


^i'h 


-0.36 - 10-3 


Ni 


-3.53 - 10-3 


Kin 


-0.19- 10-3 


^H^l-R 


-6.74- 10-3 


K. 


-0.34-10-3 


K 


-2.87- 10-3 


MLw 


0.45 - 10-3 


K 





K\.\ 


-2.14- 10-3 


Mi 


-2.31 - 10-3 


^vwN 


-0.303 


Ki 


0.13-10-3 


Mi, 


-0.38 - 10-3 






K 





Mi 


0.94 - 10-3 










Mi 


0.02 - 10-3 











Table 5: Non dimensional hydrodynamic coefficients. 
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